/*
**  File: AquaCalc.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht. 
**
**
**  AUTHOR(S) : Ton Rullmann
**  VERSION   : 1.5.0
**  DATE      : Dec 18, 1997
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"

typedef float Disterr;

typedef struct
{
    int     restr_count;
    int     num_viol;
    float   tot_viol;
    float   tot2_viol;
    float   max_viol;
}
Violsumm;

typedef struct
{
    Violsumm   *upp;
    Violsumm   *low;
    Violsumm   *lim;
}
Violrec;

typedef struct
{
    int    upp_count;
    int    low_count;
    int    uplo_count;
    int    tot_count;
}
Drstcont;


/****************  function prototypes  ****************/

#include "AquaRange.h"
#include "AquaFuncts.h"
#include "AquaStrucset.h"
#include "AquaFiles.h"
#include "MenuFuncts.h"
#include "Functs.h"

#include "AquaCalc.h"

/****************  globals  ****************/

extern Boolean  batch;                                /* import from Functs */
extern int      terseness;                            /* import from Functs */
extern Card     msg_string;                           /* import from Functs */
extern Range    *range;                               /* import from AquaRange */
char            *program_name = "AquaCalc";           /* export to Functs */
char            *program_vers = aqua_version;         /* export to Functs */

extern Filnam   dinf_filnamG, dinf_filnamNormal;      /* import from AquaFiles */

enum { opt_lim, opt_sho_selonly };

/* BEGIN_EXTH */

/****************  main  ****************/
/*
**
*/

int main( int argc, char *argv[] )
{
    char         *fct_name = "AquaCalc main";
    Strucset     theSet, selSet, *setptr;
    Model        *models = NULL;
    DistdataList inf_list, val_list;
    Distinf      *dist_inf;
    Distinfl     *dist_infl;
    Distrst      *dist_rst;
    Distval      *dist_val;
    Distvall     *dist_vall;
    DistinfTyp   ftype = inf_unktyp;
    Qtyp         qtype, qtv, qte, qti, qts;
    int          *atom_map, *residue_map;
    int          count, dist_count, tot_dist_count, a_count, s_count;
    int          struc_count, new_struc_count;
    Range        *save_range;
    Card         ident, lin;
    Filnam       dinf_filnam, drst_filnam, dval_filnam;
    Filnam       filnam_s, filnam_i, filnam_r;
    Filnam       filnam;
    FILE         *fp = NULL;
    float        lim;
    Boolean      done = false, done2 = false, err;
    int          item, how_sel, nmono;
    int          item_av = item_NULL;
    Menu         *oldmenu;
    char         temp_value[4];

    enum { item_range, item_do };
    static MenuItem menulist[] =
    {
        { item_EXEC,  "exec",  none, on },
        { item_SET,   "set",   none, on },
        { item_range, "range", none, on },
        { item_do,    "do",    none, on }, 
        { item_QUIT,  "quit",  none, on } 
    };
    static MenuItem menulist_rest_select[] =
    {
        { restsel_both,  "both",  none, on }, 
        { restsel_one,   "one",   none, on }
    };
    static MenuItem menulist_stereo[] =
    {
	{ item_average, "average",  none, on }, 
	{ item_sum,     "sum",      none, on }
    };
    static OptionItem optionlist[] =
    {
        { opt_batch,       "batch",       logical, on }, 
        { opt_terseness,   "terseness",   integer, on }, 
        { opt_prtaq,       "print",       logical, on }, 
        { opt_prtres,      "printres",    logical, on }, 
        { opt_prtmap,      "printmap",    logical, on }, 
        { opt_rest_sel,    "rest-select", oneword, on }, 
        { opt_sho_selonly, "sho-selonly", logical, on }, 
        { opt_lim,         "limit",       real,    on }
    };

    static Menu menu =
    {
        "Main menu", "command", "Give command:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    static Menu menu_rest_select =
    {
        "Rest-select menu", "type", "How to select restraints:", 
        sizeof(menulist_rest_select)/sizeof(MenuItem)
    };
    static Menu menu_stereo =
    {
        "Stereo menu", "keyword", "Give keyword:", 
        sizeof(menulist_stereo)/sizeof(MenuItem)
    };
    static Menu options =
    {
        "Options menu", "option", "Give option:", 
        sizeof(optionlist)/sizeof(MenuItem)
    }; 


/* initialize */

    ProgramInit( argc, argv );
    ProjectInit();
    setptr = &theSet;
    StrucsetInit( setptr, "theSet", true );
    SetMenu( &menu, menulist );
    SetOptions( &options, optionlist );
    SetSubMenu( &options, opt_rest_sel, &menu_rest_select, menulist_rest_select );
    ModOption( opt_lim, "0.4" );
    ModOption( opt_rest_sel, "both" );
    ModOption( opt_sho_selonly, "true" );
    sprintf( temp_value, "%d", terseness );
    ModOption( opt_terseness, temp_value );    
    range = NULL;
    inf_list.count = 0;
    val_list.count = 0;

/* the main menu loop */

    while ( true )
    {
        item = ReadMenuItem();
        switch ( item )
        {
            case item_EXEC:
		DoExec();
                break;
            case item_SET:
		ReadOption();
		batch = OptionTrue( opt_batch );
                if ( !Check_General_Options () ) {
	            printf( "Error in main: general options set are not okay\n" );
	            MError( );
                }
                terseness   = GetOptionIVal ( opt_terseness );                
                break;
            case item_range:
		oldmenu = SetOneMenu( "Give name of ranges file (or - to end):" );
		while ( !done )
		{
		    MReadWord( filnam );
		    if ( *filnam == '-' || strlen( filnam ) == 0 )
			break;
		    if ( ! ( done = OPENr( fp, filnam ) ) )
		    {
			printf( " * ranges file <%s> could not be opened\n", filnam );
			MError();
		    }
		}
		ResetMenu( oldmenu );
		if ( done ) 
		{
		    fclose( fp );
		    range = ReadRanges( filnam );
		}
                break;
            case item_do:
		goto doit;
            case item_QUIT:
            case item_NULL:
                exit( 0 );
                break;
            default:
                ProcessMenuError( menu_fatal );
                break;
        }
    }


/* do the main thing */

doit:
    
    lim = GetOptionRVal( opt_lim );

/* process the distance-values file saving the header lines and file names */

    *dval_filnam = '\0';
    ProcessDistvalFile( &val_list, &dist_count, &struc_count, &models, &qtv, 
                        setptr->filnam, dinf_filnam, drst_filnam,
		        ident, true, &ftype, dval_filnam );

    if ( ftype == inf_normal )
    {
	dist_val  = (Distval *) (val_list.lst_desc->data);
	dist_vall = NULL;
    }
    else if ( ftype == inf_extended )
    {
	dist_val  = NULL;
	dist_vall = (Distvall *) (val_list.lst_desc->data);
    }
    else {
	FATAL( "Wrong file/data type - Code error - case 3" );
        exit(99); /* Never executed but for compiler warnings */
    }

/* process the distance-restraints file (file name read from distance-values file) */

    if ( IsEmptyStr( drst_filnam ) )
	FATAL( "Name of distance-restraints file not found on distance-values file" );

    dist_rst = ProcessDistrstFile( &count, &qtype, filnam_s, filnam_i, 
                                   ident, false, drst_filnam );
    if ( NOTEQ( setptr->filnam, filnam_s ) ||
         NOTEQ( dinf_filnam, filnam_i )    || dist_count != count ) {
        FATAL( "Inconsistent information on values and restraints files" );
    }
/* process the distance-info file (file name read from distance-values file) */

    ProcessDistinfFile( &inf_list, &count, &s_count, &qtype, filnam_s, filnam_r, 
                        ident, false, &ftype, dinf_filnam );
    
    if ( ftype == inf_normal )
    {
	dist_inf  = (Distinf *) (inf_list.lst_desc->data);
	dist_infl = NULL;
    }
    else if ( ftype == inf_extended )  /* special case, cf. ProcessDistinfFile */
    {
	dist_inf  = (Distinf *) (inf_list.lst_desc->data);
	dist_infl = (Distinfl *) (inf_list.lst_desc->related->data);
    }
    else {
	FATAL( "Wrong file/data type - Code error" );
        exit(99); /* Never executed but for compiler warnings */
    }

    if ( NOTEQ( setptr->filnam, filnam_s ) ||
         NOTEQ( drst_filnam, filnam_r )    ||
         dist_count != count )
        FATAL( "Inconsistent information on values and info files" );

/* open and read the strucset file */

    if ( EMPTYSET( *setptr ) )
    {
	if ( !ProcessStrucsetFile( setptr, 0 ) )
	    exit( 1 );
    }

/* set file types */

    qte = qti = qts = undef_data; /* Against compiler warnings */
     
    switch ( qtv )
    {
	case noe_val_data:
	case noe_vall_data:
	    qte = noe_err_data;
	    qti = noe_cnt_data;
	    qts = noe_star_data;
	    break;
	case ss_val_data:
	case ss_vall_data:
	    qte = ss_err_data;
	    qti = ss_cnt_data;
	    qts = ss_star_data;
	    break;
	case hb_val_data:
	case hb_vall_data:
	    qte = hb_err_data;
	    qti = hb_cnt_data;
	    qts = hb_star_data;
	    break;
	case dis_val_data:
	case dis_vall_data:
	    qte = dis_err_data;
	    qti = dis_cnt_data;
	    qts = dis_star_data;
	    break;
	default:
	    FATAL( "Illegal qtype value for dval file - program error" );
    }


/* ----------------------------------------------- */
/* if a range was set, create a new {Distinf, Distrst, Distval} list combination */
/* use qtv to derive the restraint type */

    if ( range )
    {
	how_sel = GetOptionSubItem( opt_rest_sel );
/* create subset if necessary */
	if ( how_sel == restsel_both && OptionTrue( opt_sho_selonly ) )
	{
	    if ( ftype != inf_normal )
		FATAL( "Can't handle stereo/range combination yet; use sho-selonly off, or select restraints from file with qhow" );
	    StrucsetInit( &selSet, "selSet", true );
	    RangeStrucset( &selSet, &theSet, &atom_map, &residue_map );
	    setptr = &selSet;
	    UpdateDistInf( dist_inf, dist_count, theSet.molecule, selSet.molecule, 
	                   atom_map, residue_map );
	    WriteAquaSet( setptr, false, true );
/* warning: the subset does not contain coordinates,
 * unless ProcessStrucsetFile is called with argument -1 (see above) */
	}
        RangeDistLists( &inf_list, dist_inf, &dist_rst, &val_list, &dist_count,
	                &tot_dist_count, &a_count,  
	                &models, &new_struc_count, struc_count, setptr->molecule );
                
	if ( ftype == inf_normal ) {
	    dist_inf  = (Distinf *) (inf_list.lst_desc->data);
	    dist_infl = NULL;
	    dist_val  = (Distval *) (val_list.lst_desc->data);
	    dist_vall = NULL;
	} else if ( ftype == inf_extended ) { /* special case, cf. ProcessDistinfFile */
	    dist_inf  = (Distinf *) (inf_list.lst_desc->data);
	    dist_infl = (Distinfl *) (inf_list.lst_desc->related->data);
	    dist_val  = NULL;
	    dist_vall = (Distvall *) (val_list.lst_desc->data);
	} else {
	    FATAL( "Wrong file/data type - Code error" );
            exit(99); /* Never executed but for compiler warning */
        }
        
	struc_count = new_struc_count;
        strcpy( ident, "selected range only" );
        CreateAquaDistFiles( QtypSel( qtv ), inf_normal, "Writing files with selected restraints (normal)", true, ftype == inf_normal? true : false );
        WriteDistinfFile( dist_inf, dist_count, 1, ident, setptr->filnam );
        WriteDistrstFile( dist_rst, dist_count, ident, setptr->filnam );
	if ( ftype == inf_normal )
	{
	    WriteDistvalFile( dist_val, dist_count, ident, setptr->filnam, 
	                      models, struc_count );
	    DestroyDistdataListElem( &val_list, dist_val, true );
	}
	else if ( ftype == inf_extended )
	{
	    strcpy( dinf_filnamNormal, dinf_filnamG ); /* save the name of the original (normal) file */
	    CreateAquaDistFiles( QtypSel( qtv ), inf_extended, "Writing files with selected restraints (extended)", true, true );
	    WriteDistinflFile( dist_infl, dist_count, a_count, ident, setptr->filnam );
	    WriteDistrstFile( dist_rst, dist_count, ident, setptr->filnam );
	    WriteDistvallFile( dist_vall, dist_count, tot_dist_count, ident, 
	                       setptr->filnam, models, struc_count );
	    DestroyDistdataListElem( &val_list, dist_vall, true );
	}
	else
	    FATAL( "Wrong file/data type - Code error" );

	DestroyDistdataListElem( &inf_list, dist_inf, true );
	free( dist_rst );
	DelTxtBuffer();
	save_range = range;
	range = NULL;
	puts( "Restarting the analysis" );
	goto doit;
    }

/* ----------------------------------------------- */

    ClassifyRestraints( dist_inf, dist_count, setptr->molecule->atoms );

    switch ( ftype )
    {
	case inf_extended:
/*
	    ShowDistl( dist_infl, dist_vall, dist_rst, dist_count,
                       setptr->molecule, struc_count );
*/
	    oldmenu = SetMenu( &menu_stereo, menulist_stereo );
	    err = false;
	    while ( !done2 )
	    {
		item_av = ReadMenuItem();
		switch ( item_av )
		{
		    case item_average:
		    case item_sum:
			done2 = true;
			break;
		    case item_QUIT:
		    case item_NULL:
			err = true;
			done2 = true;
			break;
		    default:
			ProcessMenuError( menu_fatal );
			break;
		}
	    }
	    if ( item_av == item_sum )
	    {
		SetOneMenu( "Give number of monomers" );
		done2 = false;
		while ( !done2 )
		{
		    MReadWord( lin );
		    if ( sscanf( lin, "%i", &nmono ) < 1 )
		    {
			printf( "integer expected, found: %s\n", lin );
			MError();
		    }
		    else
			done2 = true;
		}
	    }
	    ResetMenu( oldmenu );
	    if ( err )
		break;
	    dist_val = DoAverage( dist_infl, dist_vall, 
	                          dist_count, struc_count, item_av, nmono );
	case inf_normal:
	    break;
	default:
	    FATAL( "Unknown dinf file type" );
    }

    DoErrors( ftype, dist_inf, dist_rst, dist_val,  
              dist_count, struc_count, setptr->molecule, models, 
	      lim, qte, qts, dval_filnam );

    DoInfo( ftype, dist_inf, dist_rst, dist_count, 
            setptr->molecule, qti, qts, dval_filnam );
    
    exit( 0 );
}

/****************  DoAverage  ****************/
/*
**  Average the distances of an extended Distval list.
**  A new, normal Distval list is created and its address returned.
**  The old list is destroyed.
**
**  NOTE: the Distval list could/should be passed back through val_list 
**  set up in main.
*/

Distval *DoAverage( Distinfl *dist_infl, Distvall *dist_vall, 
	            const int dist_count, const int struc_count, 
		    const int av_how, const int nmono )
{
    char      *fct_name = "DoAverage";
    int       d_count, s_count, count;
    Distinfl  *d_infl;
    Distvall  *d_vall;
    Distval   *d_val, *dist_val, *values;
    double    tot, minsixth = -1.0/6.0;

    if ( av_how != item_average && av_how != item_sum )
	FATAL( "Code error: av_how not average and not sum" );
    dist_val = CREATE_NEW( Distval, dist_count * struc_count );
    if ( !dist_val )
	NOMEM( "dist_val" );
    d_val   = dist_val;
    d_vall  = dist_vall;
    s_count = struc_count;

    while ( s_count-- )
    {
	d_infl  = dist_infl;
	d_count = dist_count;
	while ( d_count-- )
	{
	    values = d_vall->values;
	    count  = d_vall->count;
	    if ( count != ( d_infl->count_A * d_infl->count_B ) )
		FATAL( "inconsistent count on extended dval/dinf files" );
	    tot = 0.0;
	    while ( count-- )
	    {
		tot += pow( *values++, -6.0 );
	    }
	    if ( av_how == item_average )
		tot /= d_vall->count;
	    else
		tot /= nmono;
	    *d_val++ = pow( tot, minsixth );
	    d_infl++;
	    d_vall++;
	}

    }

    free( dist_vall->values );  /* cf. allocation scheme in ProcessDistvalFile */
    free( dist_vall );
    return( dist_val );
}

/****************  DoErrors  ****************/
/*
**  Compute and summarize the distance error values.
*/

void DoErrors( const DistinfTyp ftype, Distinf *dist_inf, Distrst *dist_rst, Distval *dist_val, 
	       const int dist_count, const int struc_count, Mol *molptr, 
	       Model *models, const float lim, Qtyp qto, Qtyp qts, Filnam dval_filnam )
{
/*    char      *fct_name = "DoErrors"; */
    FILE      *fp, *fps;
    Filnam    filnam, star_filnam;
    char      *cmmnd, *cmmnde;
    Violrec   viol_rst, viol_str, *viol_res_list, *viol_typ_list;
    Disterr   *dist_err;
/*    int       i, *typA, *typB, *typ_list; */
    int       typ_count = ContactRangeCount();


    dist_err = CalcDistError( &viol_rst, &viol_str, &viol_res_list, 
                              dist_inf, dist_rst, dist_val,  
			      dist_count, struc_count, molptr, lim );

/* open the error output file using the ident of the distance-values file */
/* write the distance error function results */

    if ( terseness < 9 ) {
        puts( "Creating the restraint violations output files" );
    }
    fp = OpenAquaFile( qto, "w", filnam, true );
    if ( !fp )
        exit( 1 );
    WriteAquaIDLines( fp, qto, filnam );
    fprintf( fp, "\n$ Distance values taken from: %s\n", dval_filnam );
    WriteTxtBuffer( fp );

    WriteDistError( &viol_rst, &viol_str, viol_res_list, 
                    dist_inf, dist_rst, dist_val, 
		    dist_err, dist_count, struc_count, 
		    molptr, models, lim, fp );

/*** second error block: test: set types to new residue number

    typA = CREATE_NEW( int, dist_count );
    typB = CREATE_NEW( int, dist_count );
    typ_list = CREATE_NEW( int, molptr->res_count );
    if ( !typA || !typB || !typ_list )
	NOMEM( "typ list(s)" );
    for ( i=0; i<dist_count; i++ )
    {
	typA[i] = dist_inf[i].res_A;
	typB[i] = dist_inf[i].res_B;
    }
    for ( i=1; i<=molptr->res_count; i++ )
	typ_list[i] = i;
    typ_count = molptr->res_count + 1;
    CalcDistErrorT2( &viol_typ_list, dist_err, typA, typB, typ_count, 
	            dist_count, struc_count, molptr, lim );
    WriteDistErrorT2( viol_typ_list, typ_list, typ_count, struc_count, 
		     models, fp );
***/

/* second error block: use types stored in dist_inf */

    CalcDistErrorT( &viol_typ_list, dist_err, dist_inf, typ_count,
	            dist_count, struc_count, molptr, lim );
    WriteDistErrorT( viol_typ_list, typ_count, struc_count, models, fp );

/* write the distance error function results to a NMR-STAR file */

    fps = OpenAquaFile( qts, "w", star_filnam, true );
    if ( !fps )
        exit( 1 );
    InitSTARFile( qts, fps );
    WriteSTARDistError( &viol_str, viol_typ_list, typ_count, struc_count, qts, fps );
    fclose( fps );

/* free the data structures */

    FreeDistError2( &viol_typ_list, struc_count );
    FreeDistError( &viol_rst, &viol_str, &viol_res_list, struc_count );
    free( dist_err );

/* clean up the error output file */

    fclose( fp );

#ifdef CLEAN0

    cmmnd = CREATE_NEW( char, 256 );
    cmmnde = GetVar( "AQUASCRIPTSDIR", true );
    sprintf( cmmnd, "%s%sclean0 %s", cmmnde, DIRSEP, filnam );
    system( cmmnd );

#endif
}

/****************  InitSTARFile  ****************/
/*
**  Initialize the NMR-STAR file with a data_label.
*/

void InitSTARFile( Qtyp qts, FILE *fp )
{
    char        *fct_name = "InitSTARFile";
    
    switch( qts )
    {
	case noe_star_data:
	    fprintf( fp, "data_AQUA_NOE\n\n" );
	    break;
	case hb_star_data:
	    fprintf( fp, "data_AQUA_HB\n\n" );
	    break;
	case ss_star_data:
	    fprintf( fp, "data_AQUA_SS\n\n" );
	    break;
	case dis_star_data:
	    fprintf( fp, "data_AQUA_DIS\n\n" );
	    break;
    default:
        FATAL("Uncoded switch");
    }    
}

/****************  DoInfo  ****************/
/*
**  Process the distance restraint information content.
*/

void DoInfo( const DistinfTyp ftype, Distinf *dist_inf, Distrst *dist_rst, 
             const int dist_count, 
             Mol *molptr, Qtyp qto, Qtyp qts, Filnam dval_filnam )
{
/*    char      *fct_name= "DoInfo"; */
    Drstcont  *cont_res, *cont_typ;
    FILE      *fp, *fps;
    Filnam    filnam, star_filnam;
    int       dtyp_count = ContactRangeCount();


    CalcDrstContent( &cont_res, &cont_typ, dist_inf, dist_rst, 
                     dist_count, dtyp_count, molptr );

/* open the output file using the ident of the distance-values file */
/* write the distance restraint information content results */

    if ( terseness < 9 ) {
        puts( "Creating the restraint statistics files" );
    }
    fp = OpenAquaFile( qto, "w", filnam, true );
    if ( !fp )
        exit( 1 );
    WriteAquaIDLines( fp, qto, filnam );
    fprintf( fp, "\n$ Distance values taken from: %s\n", dval_filnam );
    WriteTxtBuffer( fp );

    WriteDrstContent( cont_res, cont_typ, dist_count, dtyp_count, 
                      molptr, fp );

/* write the restraint statistics to an NMR-STAR file */

    fps = OpenAquaFile( qts, "a", star_filnam, true );
    if ( !fps )
        exit( 1 );
    WriteSTARDrstContent( cont_typ, dist_count, dtyp_count, qts, fps );
    fclose( fps );

}

/****************  CalcDistError  ****************/
/*
**  Calculate the distance error function.
**  A Disterr list is made, containing the violations.
**  The structure of this list is the same as of the Distval list:
**    - restraint list (Distinf/Distrst order) of structure 1
**    - restraint list of structure 2
**    - restraint list of structure 3
**    etc.
**  The violations are summarized in various Violsumm records.
*/

Disterr *CalcDistError( Violrec *viol_rst, Violrec *viol_str, Violrec **viol_res_list, 
                        Distinf *dist_inf, Distrst *dist_rst, Distval *dist_val, 
			const int dist_count, const int struc_count, Mol *molptr, 
			const float lim )
{
    char     *fct_name = "CalcDistError";
    int      d_count, s_count, idx_A, idx_B;
    Distinf  *d_inf;
    Distrst  *d_rst;
    Distval  *d_val;
    Disterr  *d_err, *e, *dist_err;
    Violsumm *v_upp, *v_low, *v_lim;
    Violrec  *viol_res, *viol_res_all;

/* compute the bound violations and store in the Disterr list */

    dist_err = CREATE_NEW( Disterr, dist_count * struc_count );
    if ( !dist_err )
	NOMEM( "dist_err" );
    d_val = dist_val;
    d_err = dist_err;
    s_count = struc_count;

    while ( s_count-- )
    {
	d_rst = dist_rst;
	d_count = dist_count;
	while ( d_count-- )
	{
	    switch ( d_rst->typ_bnd )
	    {
		case 'B':
		case 'b':
		    *d_err = *d_val - d_rst->up_bnd;     /* upper bound violation */
		    if ( *d_err < 0.F )
		    {
			*d_err = *d_val - d_rst->lo_bnd; /* lower bound violation */
							 /* (negative!) */
			if ( *d_err > 0.F )
			    *d_err = 0.F;
		    }
		    break;
		case 'u':
		    *d_err = *d_val - d_rst->up_bnd;     /* upper bound violation */
		    if ( *d_err < 0.F )
			*d_err = 0.F;
		    break;
		case 'l':
		    *d_err = *d_val - d_rst->lo_bnd;     /* lower bound violation */
						         /* (negative!) */
		    if ( *d_err > 0.F )
			*d_err = 0.F;
		    break;
		default:
		    FATAL( "program/data error: illegal bounds type" );
	    }
	    ++d_err;
	    ++d_val;
	    ++d_rst;
	}
    }

/* compute Violsumm records for all restraints */

    viol_rst->upp = NewViolsumm( dist_count );
    viol_rst->low = NewViolsumm( dist_count );
    viol_rst->lim = NewViolsumm( dist_count );

    v_upp = viol_rst->upp;
    v_low = viol_rst->low;
    v_lim = viol_rst->lim;
    d_err = dist_err;
    e     = dist_err;
    d_count = dist_count;

    while ( d_count-- )
    {
	s_count = struc_count;
	while ( s_count-- )
	{
	    UpdViolsumm( e, v_upp, v_low, v_lim, lim );
	    e += dist_count;
	}
	e = ++d_err;
	++v_upp;
	++v_low;
	++v_lim;
    }

/* compute Violsumm records for all structures and overall */

    viol_str->upp = NewViolsumm( struc_count + 1 );
    viol_str->low = NewViolsumm( struc_count + 1 );
    viol_str->lim = NewViolsumm( struc_count + 1 );

    v_upp = viol_str->upp + 1;
    v_low = viol_str->low + 1;
    v_lim = viol_str->lim + 1;
    d_err = dist_err;

    s_count = struc_count;

    while ( s_count-- )
    {
	d_count = dist_count;
	while ( d_count-- )
	{
	    UpdViolsumm( d_err, v_upp, v_low, v_lim, lim );
	    UpdViolsumm( d_err, viol_str->upp, viol_str->low, viol_str->lim, lim );
	    ++d_err;	    
	}
	++v_upp;
	++v_low;
	++v_lim;
    }

/* compute Violsumm records for all residues in array of Violrec's */
/* first Violrec: overall; remaining Violrec's: per structure */

    *viol_res_list = CREATE_NEW( Violrec, struc_count + 1 );
    if ( !viol_res_list )
	NOMEM( "viol_res_list" );
    viol_res_all = *viol_res_list;
    viol_res_all->upp = NewViolsumm( molptr->res_count );
    viol_res_all->low = NewViolsumm( molptr->res_count );
    viol_res_all->lim = NewViolsumm( molptr->res_count );

    viol_res = viol_res_all + 1;
    s_count = struc_count;
    while ( s_count-- )
    {
	viol_res->upp = NewViolsumm( molptr->res_count );
	viol_res->low = NewViolsumm( molptr->res_count );
	viol_res->lim = NewViolsumm( molptr->res_count );
	viol_res++;
    }
    
    viol_res = viol_res_all + 1;
    d_err = dist_err;

    s_count = struc_count;

    while ( s_count-- )
    {
	d_count = dist_count;
	d_inf = dist_inf;
	while ( d_count-- )
	{
	    idx_A = FindResidueIndexFromNum( d_inf->res_A, molptr );
	    idx_B = FindResidueIndexFromNum( d_inf->res_B, molptr );
	    if ( idx_A < 0 || idx_B < 0 )
		FATAL( "Indices out of range - Code error" );
	    UpdViolsumm( d_err, viol_res->upp+idx_A, viol_res->low+idx_A, 
	                        viol_res->lim+idx_A, lim );
	    UpdViolsumm( d_err, viol_res_all->upp+idx_A, viol_res_all->low+idx_A, 
	                        viol_res_all->lim+idx_A, lim );
	    if ( idx_B != idx_A )
	    {
		UpdViolsumm( d_err, viol_res->upp+idx_B, viol_res->low+idx_B, 
				    viol_res->lim+idx_B, lim );
		UpdViolsumm( d_err, viol_res_all->upp+idx_B, viol_res_all->low+idx_B, 
				    viol_res_all->lim+idx_B, lim );
	    }
	    ++d_err;
	    ++d_inf;
	}
	++viol_res;
    }

    return( dist_err );
}

/****************  CalcDistErrorT  ****************/
/*
**  Calculate the distance error function (continued).
**  A Disterr list, containing the violations, should have been made
**  by CalcDistError.
**  The violations are here summarized in Violsumm records, categorizing
**  on types (cf. categorization by residue in CalcDistError).
**  For the types it is assumed that dist_inf->dtype has been set to a
**  ContactRange value, which should be between 0 and typ_count-1.
*/

void CalcDistErrorT( Violrec **viol_typ_list, Disterr *dist_err, 
                     Distinf *dist_inf, const int typ_count, 
	             const int dist_count, const int struc_count, Mol *molptr, 
		     const float lim )
{
    char     *fct_name = "CalcDistErrorT";
    int      d_count, s_count, typ;
    Distinf  *d_inf;
    Disterr  *d_err;
    Violrec  *viol_typ, *viol_typ_all;

/* compute Violsumm records for all types in array of Violrec's */
/* first Violrec: overall; remaining Violrec's: per structure */

    *viol_typ_list = CREATE_NEW( Violrec, struc_count + 1 );
    if ( !viol_typ_list )
	NOMEM( "viol_typ_list" );
    viol_typ_all = *viol_typ_list;
    viol_typ_all->upp = NewViolsumm( typ_count );
    viol_typ_all->low = NewViolsumm( typ_count );
    viol_typ_all->lim = NewViolsumm( typ_count );

    viol_typ = viol_typ_all + 1;
    s_count = struc_count;
    while ( s_count-- )
    {
	viol_typ->upp = NewViolsumm( typ_count );
	viol_typ->low = NewViolsumm( typ_count );
	viol_typ->lim = NewViolsumm( typ_count );
	viol_typ++;
    }
    
    viol_typ = viol_typ_all + 1;
    d_err = dist_err;

    s_count = struc_count;

    while ( s_count-- )
    {
	d_count = dist_count;
	d_inf = dist_inf;
	while ( d_count-- )
	{
	    typ = d_inf->dtype;
	    if ( typ < 0 || typ >= typ_count )
	    {
		printf( "%i %i\n", typ, typ_count );
		FATAL( "Type index out of range - Code error" );
	    }
	    UpdViolsumm( d_err, viol_typ->upp+typ, viol_typ->low+typ, 
	                        viol_typ->lim+typ, lim );
	    UpdViolsumm( d_err, viol_typ_all->upp+typ, viol_typ_all->low+typ, 
	                        viol_typ_all->lim+typ, lim );
	    ++d_err;
	    ++d_inf;
	}
	++viol_typ;
    }
}

/****************  CalcDistErrorT2  ****************/
/*
**  Calculate the distance error function (continued).
**  A Disterr list, containing the violations, should have been made
**  by CalcDistError.
**  The violations are here summarized in Violsumm records, categorizing
**  on types (cf. categorization by residue in CalcDistError).
**  In this version we have an atom-type (two indices, stored in typA and typB).
*/

void CalcDistErrorT2( Violrec **viol_typ_list, Disterr *dist_err, 
                     int *typA, int *typB, const int typ_count, 
	             const int dist_count, const int struc_count, Mol *molptr, 
		     const float lim )
{
    char     *fct_name = "CalcDistErrorT2";
    int      d_count, s_count, idx_A, idx_B, *tA, *tB;
    Disterr  *d_err;
    Violrec  *viol_typ, *viol_typ_all;

/* compute Violsumm records for all types in array of Violrec's */
/* first Violrec: overall; remaining Violrec's: per structure */

    *viol_typ_list = CREATE_NEW( Violrec, struc_count + 1 );
    if ( !viol_typ_list )
	NOMEM( "viol_typ_list" );
    viol_typ_all = *viol_typ_list;
    viol_typ_all->upp = NewViolsumm( typ_count );
    viol_typ_all->low = NewViolsumm( typ_count );
    viol_typ_all->lim = NewViolsumm( typ_count );

    viol_typ = viol_typ_all + 1;
    s_count = struc_count;
    while ( s_count-- )
    {
	viol_typ->upp = NewViolsumm( typ_count );
	viol_typ->low = NewViolsumm( typ_count );
	viol_typ->lim = NewViolsumm( typ_count );
	viol_typ++;
    }
    
    viol_typ = viol_typ_all + 1;
    d_err = dist_err;

    s_count = struc_count;

    while ( s_count-- )
    {
	d_count = dist_count;
	tA = typA;
	tB = typB;
	while ( d_count-- )
	{
	    idx_A = *tA++;
	    idx_B = *tB++;
	    if ( idx_A < 0 || idx_B < 0 || idx_A >= typ_count || idx_B >= typ_count )
	    {
		printf( "%i %i %i\n", idx_A, idx_B, typ_count );
		FATAL( "Indices out of range - Code error" );
	    }
	    UpdViolsumm( d_err, viol_typ->upp+idx_A, viol_typ->low+idx_A, 
	                        viol_typ->lim+idx_A, lim );
	    UpdViolsumm( d_err, viol_typ_all->upp+idx_A, viol_typ_all->low+idx_A, 
	                        viol_typ_all->lim+idx_A, lim );
	    if ( idx_B != idx_A )
	    {
		UpdViolsumm( d_err, viol_typ->upp+idx_B, viol_typ->low+idx_B, 
				    viol_typ->lim+idx_B, lim );
		UpdViolsumm( d_err, viol_typ_all->upp+idx_B, viol_typ_all->low+idx_B, 
				    viol_typ_all->lim+idx_B, lim );
	    }
	    ++d_err;
	}
	++viol_typ;
    }
}

/****************  WriteDistError  ****************/
/*
**  Write the results of the distance error calculation.
*/

void WriteDistError( Violrec *viol_rst, Violrec *viol_str, Violrec *viol_res, 
                     Distinf *dist_inf, Distrst *dist_rst, Distval *dist_val, 
		     Disterr *dist_err, const int dist_count, const int struc_count,
		     Mol *molptr, Model *models, const float lim, FILE *fp)
{
    int      d_count, s_count, r_count, res_count;
    Distinf  *d_inf;
    Distrst  *d_rst;
    Distval  *d_val, *v;
    Disterr  *d_err, *e;
    Violsumm *v_upp, *v_low, *v_lim;
    Atom     *atoms;
    Residue  *residues;
    Model    *mdl;

    atoms = molptr->atoms;
    res_count = molptr->res_count;
    
    fputc( '\n', fp );
    fprintf( fp, "$ EXPLAIN\n" );
    fprintf( fp, "$ NumUPPv: number of upper bound violations\n" );
    fprintf( fp, "$ MaxUPPv: largest upper bound violation\n" );
    fprintf( fp, "$ TotUPPv: total value of upper bound violations\n" );
    fprintf( fp, "$ Av1UPPv: average value of upper bound violations (calc. over violated bounds)\n" );
    fprintf( fp, "$ sd1UPPv: standard deviation of upper bound violations (calc. over violated bounds)\n" );
    fprintf( fp, "$ Av2UPPv: average value of upper bound violations (calc. over all restraints)\n" );
    fprintf( fp, "$ sd2UPPv: standard deviation of upper bound violations (calc. over all restraints)\n" );
    fprintf( fp, "$ ...LOWv: idem for lower bound violations\n" );
    fprintf( fp, "$ ...BIGv: idem for big upper and lower bound violations (i.e. abs. values > %5.2f)\n", lim );
    fprintf( fp, "$ NRestr : number of restraints used in summation\n" );
    fprintf( fp, "$ END\n" );

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  GENERAL\n" );
    fprintf( fp, "$ LIST #RECORDS %i\n", 3 );
    fprintf( fp, "NRestraints %i\n", dist_count );
    fprintf( fp, "NStructures %i\n", struc_count );
    fprintf( fp, "NResidues   %i\n", res_count );
    fprintf( fp, "$ END\n" );

/* write the distances */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  DISTANCES\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", dist_count, 11 + struc_count );
    fprintf( fp, "Chn Resnm Res Atom  Chn Resnm Res Atom  Type  LoBnd UpBnd" );
    s_count = struc_count;
    mdl = models;
    while ( s_count-- )
	fprintf( fp, " Str_%02i", mdl++->struc_num );
    fputc( '\n', fp );

    d_inf = dist_inf;
    d_rst = dist_rst;
    d_val = dist_val;
    v     = dist_val;
    d_count = dist_count;
    while ( d_count-- )
    {
	WriteFAtom( fp, &atoms[ d_inf->idx_A ], use_ext );
	fputc( ' ', fp );
	WriteFAtom( fp, &atoms[ d_inf->idx_B ], use_ext );
	fprintf( fp, " %3i %6.2f %6.2f",
		d_inf->dtype, d_rst->lo_bnd, d_rst->up_bnd );
	s_count = struc_count;
	while ( s_count-- )
	{
	    fprintf( fp, " %6.2f", *v );
	    v += dist_count;
	}
	fputc( '\n', fp );
	v = ++d_val;
	++d_inf;
	++d_rst;
    }
    fprintf( fp, "$ END\n" );

/* write the restraint violations */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATIONS\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", dist_count, 11 + struc_count );
    fprintf( fp, "Chn Resnm Res Atom  Chn Resnm Res Atom  Type  LoBnd UpBnd" );
    s_count = struc_count;
    mdl = models;
    while ( s_count-- )
	fprintf( fp, " Str_%02i", mdl++->struc_num );
    fputc( '\n', fp );

    d_inf = dist_inf;
    d_rst = dist_rst;
    d_err = dist_err;
    e     = dist_err;
    d_count = dist_count;
    while ( d_count-- )
    {
	WriteFAtom( fp, &atoms[ d_inf->idx_A ], use_ext );
	fputc( ' ', fp );
	WriteFAtom( fp, &atoms[ d_inf->idx_B ], use_ext );
	fprintf( fp, " %3i %6.2f %6.2f",
		d_inf->dtype, d_rst->lo_bnd, d_rst->up_bnd );
	s_count = struc_count;
	while ( s_count-- )
	{
	    fprintf( fp, " %6.2f", *e );
	    e += dist_count;
	}
	fputc( '\n', fp );
	e = ++d_err;
	++d_inf;
	++d_rst;
    }
    fprintf( fp, "$ END\n" );

/* write the overall restraint violation summary */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_OVERALL\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", 1, 1 + 18 );
    fprintf( fp, "Structure" );
    WriteHdrViolsummAv( fp );
    v_upp = viol_str->upp;
    v_low = viol_str->low;
    v_lim = viol_str->lim;
    fprintf( fp, "Overall  " );
    WriteViolsummAv( fp, v_upp, v_low, v_lim );
    fprintf( fp, "$ END\n" );

/* write the restraint violation summary per structure */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_STRUCTURE\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", struc_count, 1 + 18 );
    fprintf( fp, "Structure" );
    WriteHdrViolsummAv( fp );
    v_upp = viol_str->upp + 1;
    v_low = viol_str->low + 1;
    v_lim = viol_str->lim + 1;
    s_count = struc_count;
    mdl = models;
    while ( s_count-- )
    {
	fprintf( fp, "Str_%02i   ", mdl++->struc_num );
	WriteViolsummAv( fp, v_upp++, v_low++, v_lim++ );
    }
    fprintf( fp, "$ END\n" );

/* write the restraint violation summary per restraint */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_RESTRAINT\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", dist_count, 8 + 10 );
    fprintf( fp, "Chn Resnm Res Atom  Chn Resnm Res Atom " );
    WriteHdrViolsumm( fp );

    d_inf = dist_inf;
    d_rst = dist_rst;
    v_upp = viol_rst->upp;
    v_low = viol_rst->low;
    v_lim = viol_rst->lim;
    d_count = dist_count;
    while ( d_count-- )
    {
	WriteFAtom( fp, &atoms[ d_inf->idx_A ], use_ext );
	fputc( ' ', fp );
	WriteFAtom( fp, &atoms[ d_inf->idx_B ], use_ext );
	WriteViolsumm( fp, v_upp++, v_low++, v_lim++ );
	++d_inf;
	++d_rst;
    }
    fprintf( fp, "$ END\n" );

/* write the restraint violation summary per residue */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_RESIDUE_for_ALL_STRUCTURES\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", res_count, 3 + 10 );
    fprintf( fp, "Chn Resnm Res " );
    WriteHdrViolsumm( fp );

    v_upp = viol_res->upp;
    v_low = viol_res->low;
    v_lim = viol_res->lim;
    r_count = 1;
    residues = molptr->residues;
    while ( r_count++ <= res_count )
    {
	WriteFResidue( fp, residues++ );
	WriteViolsumm( fp, v_upp++, v_low++, v_lim++ );
    }
    fprintf( fp, "$ END\n" );

    viol_res++;
    s_count = struc_count;
    mdl = models;
    while ( s_count-- )
    {
	fputc( '\n', fp );
	fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_RESIDUE_for_STRUCTURE_Str_%02i\n", mdl++->struc_num );
	fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", res_count, 3 + 10 );
	fprintf( fp, "Chn Resnm Res " );
	WriteHdrViolsumm( fp );
    
	v_upp = viol_res->upp;
	v_low = viol_res->low;
	v_lim = viol_res->lim;
	r_count = 1;
	residues = molptr->residues;
	while ( r_count++ <= res_count )
	{
	    WriteFResidue( fp, residues++ );
	    WriteViolsumm( fp, v_upp++, v_low++, v_lim++ );
	}
	fprintf( fp, "$ END\n" );
	viol_res++;
    }
}

/****************  WriteDistErrorT  ****************/
/*
**  Write the results of the distance error calculation per type.
**  For the types it is assumed that they range between 0 and typ_count-1, 
**  (skipping the first, i.e. the undefined category in the output)
**  and the type names are delivered by ContactRangeNam (cf. CalcDistErrorT).
*/

void WriteDistErrorT( Violrec *viol_typ, const int typ_count, const int struc_count, 
		     Model *models, FILE *fp )
{
    int      s_count, t_count;
    Violsumm *v_upp, *v_low, *v_lim;
    Model    *mdl;

/* write the restraint violation summary per type */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_RANGE_TYPE_for_ALL_STRUCTURES\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", typ_count-1, 1 + 18 );
    fprintf( fp, "Type       " );
    WriteHdrViolsummAv( fp );

    v_upp = viol_typ->upp+1;
    v_low = viol_typ->low+1;
    v_lim = viol_typ->lim+1;
    t_count = 1;
    while ( t_count < typ_count )
    {
	fprintf( fp, "%10s ", ContactRangeNam( (ContactRange) (t_count++) ) );
	WriteViolsummAv( fp, v_upp++, v_low++, v_lim++ );
    }
    fprintf( fp, "$ END\n" );

    viol_typ++;
    s_count = struc_count;
    mdl = models;
    while ( s_count-- )
    {
	fputc( '\n', fp );
	fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_RANGE_TYPE_for_STRUCTURE_Str_%02i\n", mdl++->struc_num );
	fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", typ_count-1, 1 + 18 );
	fprintf( fp, "Type       " );
	WriteHdrViolsummAv( fp );
    
	v_upp = viol_typ->upp+1;
	v_low = viol_typ->low+1;
	v_lim = viol_typ->lim+1;
	t_count = 1;
	while ( t_count < typ_count )
	{
	    fprintf( fp, "%10s ", ContactRangeNam( (ContactRange) (t_count++) ) );
	    WriteViolsummAv( fp, v_upp++, v_low++, v_lim++ );
	}
	fprintf( fp, "$ END\n" );
	viol_typ++;
    }
}

/****************  WriteDistErrorT2  ****************/
/*
**  Write the results of the distance error calculation per type.
*/

void WriteDistErrorT2( Violrec *viol_typ, int *typ_list, 
                     const int typ_count, const int struc_count, 
		     Model *models, FILE *fp )
{
    int      s_count, t_count, *typ;
    Violsumm *v_upp, *v_low, *v_lim;
    Model    *mdl;

/* write the restraint violation summary per type */

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_TYPE_for_ALL_STRUCTURES\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", typ_count, 1 + 10 );
    fprintf( fp, "Type " );
    WriteHdrViolsumm( fp );

    v_upp = viol_typ->upp;
    v_low = viol_typ->low;
    v_lim = viol_typ->lim;
    t_count = 1;
    typ = typ_list;
    while ( t_count++ <= typ_count )
    {
	fprintf( fp, "%4i ", *typ++ );
	WriteViolsumm( fp, v_upp++, v_low++, v_lim++ );
    }
    fprintf( fp, "$ END\n" );

    viol_typ++;
    s_count = struc_count;
    mdl = models;
    while ( s_count-- )
    {
	fputc( '\n', fp );
	fprintf( fp, "$ DATA  RESTRAINT_VIOLATION_SUMMARY_per_TYPE_for_STRUCTURE_Str_%02i\n", mdl++->struc_num );
	fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", typ_count, 1 + 10 );
	fprintf( fp, "Type " );
	WriteHdrViolsumm( fp );
    
	v_upp = viol_typ->upp;
	v_low = viol_typ->low;
	v_lim = viol_typ->lim;
	t_count = 1;
	typ = typ_list;
	while ( t_count++ <= typ_count )
	{
	    fprintf( fp, "%4i ", *typ++ );
	    WriteViolsumm( fp, v_upp++, v_low++, v_lim++ );
	}
	fprintf( fp, "$ END\n" );
	viol_typ++;
    }
}


/****************  WriteSTARDistError  ****************/
/*
**  Write the results of the distance error calculation in NMR-STAR format.
**  Updated NMR-STAR tags and values, written by JFD @ Tue May  9 12:16:42 CDT 2000 
*/                 

void WriteSTARDistError( Violrec *viol_str, Violrec *viol_typ, const int typ_count, 
                         const int struc_count, Qtyp qts, FILE *fp )
#define size_KEY 20
{
    char      *fct_name = "WriteSTARDistError";
    Violsumm  *v_upp, *v_low;
    float     av_upp, rms_upp, av_low, rms_low;
    char      *key;
/*    int       t_count;
    char      idtxt[size_KEY]; */

    char   *Mol_label = GetVar( "AQUAPROJECT", false );

    switch( qts )
    {
	case noe_star_data:
	    key = "NOE";
	    break;
	case hb_star_data:
	    key = "HB";
	    break;
	case ss_star_data:
	    key = "SS";
	    break;
	case dis_star_data:
	    key = "DIS";
	    break;
	default:
            FATAL("Unknow case; code error ?");
	    exit(99); /* Will never be executed */
    }

    if ( qts != noe_star_data && qts != hb_star_data )
	return;    

    fprintf( fp, "save_%s_constraint_violation_stats\n", key );
    fprintf( fp, "   _Saveframe_category   distance_constraint_violation_stats\n\n" );

    fprintf( fp, "   _Details\n;\n" );
    fprintf( fp, "Violation averages and rmsd are calculated over all restraints and all models.\n" );
    fprintf( fp, ";\n" );

    fprintf( fp, "   loop_\n" );
    fprintf( fp, "       _Constraint_list_label\n" );
    fprintf( fp, "       $distance_constraint_%s_%s\n", key, Mol_label );
    fprintf( fp, "\n   stop_\n\n" );
    
    fprintf( fp, "   _Conformer_total_number                    %3d\n\n", struc_count );

    /* write the overall restraint violation summary */

    if ( qts == noe_star_data )
    {    
	v_upp = viol_str->upp;
	v_low = viol_str->low;

/* Not used now, check before reusing...
	v_lim = viol_str->lim;
        
	CombViolsumm( v_upp, v_low, &v_comb );
*/
              
/* JFD @ Wed Jul 26 11:01:41 CDT 2000
   added another rms calculation
*/
        RmsViolsumm2( v_upp,  &av_upp, &rms_upp, false );
        RmsViolsumm2( v_low,  &av_low, &rms_low, false );

/*	RmsViolsumm2( &v_comb, &avB,    &rmsB,    false ); */
        
/*    Not using the combination of the two for now.
	fprintf( fp, "  _Total_distance_constraint_violation_max " );
	fprintf( fp, " %7.2f\n", v_comb.max_viol );
	fprintf( fp, "  _Total_distance_constraint_violation_avg " );
	fprintf( fp, " %7.3f\n", avB );
	fprintf( fp, "  _Total_distance_constraint_violation_rmsd" );
	fprintf( fp, " %7.3f\n", rmsB );
*/
        
	fprintf( fp, "   _NOE_upper_violation_max " );
	fprintf( fp, " %9.4f\n", v_upp->max_viol );
	fprintf( fp, "   _NOE_upper_violation_rms " );
	fprintf( fp, " %9.4f\n", rms_upp );
	fprintf( fp, "   _NOE_upper_violation_av  " );
	fprintf( fp, " %9.4f\n", av_upp );

        
	fprintf( fp, "   _NOE_lower_violation_max " );
	fprintf( fp, " %9.4f\n", -(v_low->max_viol) );
	fprintf( fp, "   _NOE_lower_violation_rms " );
	fprintf( fp, " %9.4f\n", rms_low );
	fprintf( fp, "   _NOE_lower_violation_av  " );
	fprintf( fp, " %9.4f\n", -(av_low) );

    }

    fprintf( fp, "\nsave_\n\n" );
    
/* some old code for this routine:   
    if ( qts == noe_star_data )
    {
	v_upp = viol_typ->upp+1;
	v_low = viol_typ->low+1;
	v_lim = viol_typ->lim+1;
	t_count = 1;
	while ( t_count < typ_count )
	{
	    strcpy( idtxt, "'" );
	    strcat( idtxt, ContactRangeStarNam( t_count++ ) );
	    strcat( idtxt, "'" );
	    fprintf( fp, "%-30s", idtxt );
	    CombViolsumm( v_upp++, v_low++, &v_comb );
	    AvRmsViolsumm( &v_comb, &avB, &rmsB, false );
	    fprintf( fp, " %7.3f %7.3f    .\n", avB, rmsB );
	}
    }
    else if ( qts == hb_star_data )
    {
	v_upp = viol_str->upp;
	v_low = viol_str->low;
	v_lim = viol_str->lim;
    
	CombViolsumm( v_upp, v_low, &v_comb );
	AvRmsViolsumm( &v_comb, &avB, &rmsB, false );
	
	fprintf( fp, "%-30s %7.3f %7.3f    .\n", "'hydrogen bond'", avB, rmsB );
    }
*/
        
    
}
#undef size_KEY

/****************  FreeDistError  ****************/
/*
**  Free the memory allocated in CalcDistError.
*/

void FreeDistError( Violrec *viol_rst, Violrec *viol_str, 
                    Violrec **viol_res_list, const int struc_count )
{
    FreeDistError2( viol_res_list, struc_count );

    free( viol_str->upp );
    free( viol_str->low );
    free( viol_str->lim );

    free( viol_rst->upp );
    free( viol_rst->low );
    free( viol_rst->lim );
}

/****************  FreeDistError2  ****************/
/*
**  Free the memory allocated in CalcDistError[2] for
**  viol_res_list or viol_typ_list.
*/

void FreeDistError2( Violrec **viol_x_list, const int struc_count )
{
    Violrec  *viol_x;
    int      count;

    viol_x = *viol_x_list;
    count = struc_count + 1;
    while ( count-- )
    {
	free( viol_x->upp );
	free( viol_x->low );
	free( viol_x->lim );
	viol_x++;
    }
    
    free( *viol_x_list );
}

/****************  NewViolsumm  ****************/
/*
**  Allocate and initialize a Violsumm array.
*/

Violsumm *NewViolsumm( int count )
{
    Violsumm   *new_variable, *new_array;

    new_array = CREATE_NEW( Violsumm, count );
    if ( !new_array )
	NoMemErr( "NewViolsumm" );

    new_variable = new_array;
    while ( count-- )
    {
	new_variable->restr_count = 0;
	new_variable->num_viol = 0;
	new_variable->tot_viol = 0.F;
	new_variable->tot2_viol = 0.F;
	new_variable->max_viol = 0.F;
	new_variable++;
    }
    return( new_array );
}

/****************  UpdViolsumm  ****************/
/*
**  Update a set of Violsumm structs (v_upp, v_low, v_lim) using the error
**  value *e.
*/

void UpdViolsumm( float *e, Violsumm *v_upp, Violsumm *v_low, Violsumm *v_lim, 
                  const float lim )
{
    float   abse;

    if ( *e > 0.F )
    {
	abse = *e;
	v_upp->num_viol += 1;
	v_upp->tot_viol += *e;
	v_upp->tot2_viol += ( abse * abse );
	if ( *e > v_upp->max_viol )
	    v_upp->max_viol = *e;
    }
    else if ( *e < 0.F )
    {
	abse = -(*e);
	v_low->num_viol += 1;
	v_low->tot_viol += *e;
	v_low->tot2_viol += ( abse * abse );
	if ( *e < v_low->max_viol )
	    v_low->max_viol = *e;
    }
    else
	abse = 0.F;

    if ( abse > lim )
    {
	v_lim->num_viol += 1;
	v_lim->tot_viol += abse;
	v_lim->tot2_viol += ( abse * abse );
	if ( abse > v_lim->max_viol )
	    v_lim->max_viol = abse;
    }

    v_upp->restr_count += 1;
    v_low->restr_count += 1;
    v_lim->restr_count += 1;
}

/****************  AvRmsViolsumm  ****************/
/*
**  Compute average and rmsd for a Violsumm struct.
*/

void  AvRmsViolsumm( Violsumm *v, float *av, float *rms, const Boolean numViol )
{
    int   n;
    
    if ( numViol )
        n = v->num_viol;
    else
        n = v->restr_count;
    if ( n == 0 )
    {
	*av = 0.F;
	*rms = 0.F;
    }
    else
    {
	*av = v->tot_viol / n;
        *rms = sqrt ( v->tot2_viol / n - ( *av * *av ) );
    }
}


/****************  RmsViolsumm2  ****************/
/*
**  Compute different rmsd for a Violsumm struct.
*/

void  RmsViolsumm2( Violsumm *v, float *av, float *rms, const Boolean numViol )
{
    int   n;
    
    if ( numViol )
        n = v->num_viol;
    else
        n = v->restr_count;
    if ( n == 0 )
    {
	*av = 0.F;
	*rms = 0.F;
    }
    else
    {
	*av = v->tot_viol / n;
        *rms = sqrt ( v->tot2_viol / n );
    }
}


/****************  CombViolsumm  ****************/
/*
**  Combine the upper and lower bound data for a Violsumm struct.
*/

void CombViolsumm( Violsumm *v_upp, Violsumm *v_low, Violsumm *v_comb )
{
    char  *fct_name = "CombViolsumm";

    if ( v_upp->restr_count != v_low->restr_count )
	FATAL( "unequal restraint counts - code error" );
    v_comb->restr_count = v_upp->restr_count;
    v_comb->num_viol    = v_upp->num_viol    + v_low->num_viol;  /* CORRECT? */
    v_comb->tot_viol    = v_upp->tot_viol    - v_low->tot_viol;
    v_comb->tot2_viol   = v_upp->tot2_viol   + v_low->tot2_viol;
    v_comb->max_viol    = MAX( v_upp->max_viol, -(v_low->max_viol) );
}

/****************  WriteViolsumm  ****************/
/*
**  Write a set of Violsumm structs (v_upp, v_low, v_lim).
*/

void WriteViolsumm( FILE *fp, Violsumm *v_upp, Violsumm *v_low, Violsumm *v_lim )
{
    fprintf( fp, " %7i %7.2f %7.2f %7i %7.2f %7.2f %7i %7.2f %7.2f  %6i\n", 
	     v_upp->num_viol, v_upp->max_viol, v_upp->tot_viol, 
	     v_low->num_viol, v_low->max_viol, v_low->tot_viol, 
	     v_lim->num_viol, v_lim->max_viol, v_lim->tot_viol, 
	     v_lim->restr_count );
}

/****************  WriteViolsummAv  ****************/
/*
**  Write a set of Violsumm structs (v_upp, v_low, v_lim), 
**  including averages info.
*/

void WriteViolsummAv( FILE *fp, Violsumm *v_upp, Violsumm *v_low, Violsumm *v_lim )
{
    float   avA, rmsA, avB, rmsB;
    
    AvRmsViolsumm( v_upp, &avA, &rmsA, true );
    AvRmsViolsumm( v_upp, &avB, &rmsB, false );
    fprintf( fp, " %7i %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f", 
	     v_upp->num_viol, v_upp->max_viol, v_upp->tot_viol, 
	     avA, rmsA, avB, rmsB );
    AvRmsViolsumm( v_low, &avA, &rmsA, true );
    AvRmsViolsumm( v_low, &avB, &rmsB, false );
    fprintf( fp, " %7i %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f", 
	     v_low->num_viol, v_low->max_viol, v_low->tot_viol, 
	     avA, rmsA, avB, rmsB );
    fprintf( fp, " %7i %7.2f %7.2f  %6i\n", 
	     v_lim->num_viol, v_lim->max_viol, v_lim->tot_viol, 
	     v_lim->restr_count );
}

/****************  WriteHdrViolsumm  ****************/
/*
**  Write a header for subsequent calls to WriteViolsumm.
*/

void WriteHdrViolsumm( FILE *fp )
{
    fprintf( fp, 
    " NumUPPv MaxUPPv TotUPPv NumLOWv MaxLOWv TotLOWv NumBIGv MaxBIGv TotBIGv  NRestr\n" );
}

/****************  WriteHdrViolsummAv  ****************/
/*
**  Write a header for subsequent calls to WriteViolsummAv.
*/

void WriteHdrViolsummAv( FILE *fp )
{
    fprintf( fp, 
    " NumUPPv MaxUPPv TotUPPv Av1UPPv sd1UPPv Av2UPPv sd2UPPv" );
    fprintf( fp, 
    " NumLOWv MaxLOWv TotLOWv Av1LOWv sd1LOWv Av2LOWv sd2LOWv" );
    fprintf( fp, 
    " NumBIGv MaxBIGv TotBIGv  NRestr\n" );
}

/****************  CalcDrstContent  ****************/
/*
**  Calculate the distance restraint information content.
**  For the types it is assumed that dist_inf->dtype has been set to a
**  ContactRange value, which should be between 0 and dtyp_count-1.
*/

void CalcDrstContent( Drstcont **cont_res, Drstcont **cont_typ, 
                      Distinf *dist_inf, Distrst *dist_rst, 
                      const int dist_count, const int dtyp_count, Mol *molptr )
{
    char     *fct_name = "CalcDrstContent";
    int      d_count, typ, idx_A, idx_B;
    Distinf  *d_inf;
    Distrst  *d_rst;

/* calculate the number of constraints per residue and per type */

    *cont_res = NewDrstcont( molptr->res_count );
    *cont_typ = NewDrstcont( dtyp_count );

    d_count = dist_count;
    d_inf = dist_inf;
    d_rst = dist_rst;
    while ( d_count-- )
    {
	idx_A = FindResidueIndexFromNum( d_inf->res_A, molptr );
	idx_B = FindResidueIndexFromNum( d_inf->res_B, molptr );
	if ( idx_A < 0 || idx_B < 0 )
	    FATAL( "Indices out of range - Code error" );
	UpdDrstcont( *cont_res + idx_A, d_rst );
	if ( idx_B != idx_A )
	    UpdDrstcont( *cont_res + idx_B, d_rst );

	typ = d_inf->dtype;
	if ( typ < 0 || typ >= dtyp_count )
	    FATAL( "Type index out of range - Code error" );
	UpdDrstcont( *cont_typ + typ, d_rst );
	++d_inf;
	++d_rst;
    }
}

/****************  WriteDrstContent  ****************/
/*
**  Write the distance restraint information content.
**  For the types it is assumed that they range between 0 and dtyp_count-1, 
**  (skipping the first, i.e. the undefined category in the output)
**  and the type names are delivered by ContactRangeNam (cf. CalcDrstContent).
*/

void WriteDrstContent( Drstcont *cont_res, Drstcont *cont_typ,
                       const int dist_count, const int dtyp_count, 
		       Mol *molptr, FILE *fp )
{
    int      t_count, r_count, res_count;
    Residue  *residues;
    
    res_count = molptr->res_count;

    fputc( '\n', fp );   
    fprintf( fp, "$ EXPLAIN\n" );
    fprintf( fp, "$ NLower: number of restraints with only lower bound\n" );
    fprintf( fp, "$ NUpper: number of restraints with only upper bound\n" );
    fprintf( fp, "$ NUp&Lo: number of restraints with upper and lower bound\n" );
    fprintf( fp, "$ NTotal: total number of restraints\n" );
    fprintf( fp, "$ END\n" );

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  GENERAL\n" );
    fprintf( fp, "$ LIST #RECORDS %i\n", 2 );
    fprintf( fp, "NRestraints %i\n", dist_count );
    fprintf( fp, "NResidues   %i\n", res_count );
    fprintf( fp, "$ END\n" );
    
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  NUMBER_OF_RESTRAINTS_per_RESIDUE\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", res_count, 3 + 4 );
    fprintf( fp, "Chn Resnm Res " );
    WriteHdrDrstcont( fp );

    r_count = 1;
    residues = molptr->residues;
    while ( r_count++ <= res_count )
    {
	WriteFResidue( fp, residues++ );
	WriteDrstcont( fp, cont_res++ );
    }
    fprintf( fp, "$ END\n" );

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  NUMBER_OF_RESTRAINTS_per_RANGE_TYPE\n" );
    fprintf( fp, "$ TABLE #RECORDS %i #ITEMS %i\n", dtyp_count-1, 1 + 4 );
    fprintf( fp, "Type       " );
    WriteHdrDrstcont( fp );

    t_count = 1;
    cont_typ++;
    while ( t_count < dtyp_count )
    {
	fprintf( fp, "%10s ", ContactRangeNam( (ContactRange) (t_count++) ) );
	WriteDrstcont( fp, cont_typ++ );
    }
    fprintf( fp, "$ END\n" );
}

/****************  WriteSTARDrstContent  ****************/
/*
**  Write the distance restraint information content in NMR-STAR format.
*/

void WriteSTARDrstContent( Drstcont *cont_typ,
                           const int dist_count, const int dtyp_count, 
		           Qtyp qts, FILE *fp )
{
    char   *fct_name = "WriteSTARDrstContent";
    char   *key, *Mol_label = GetVar( "AQUAPROJECT", false );
    int    ix;

    switch( qts )
    {
	case noe_star_data:
	    key = "NOE";
	    fprintf( fp, "save_NOE_constraint_number_stats\n" );
	    break;
	case hb_star_data:
	    key = "HB";
	    fprintf( fp, "save_hydrogen_bond_number_stats\n" );
	    break;
	case ss_star_data:
	    key = "SS";
	    fprintf( fp, "save_disulfide_bond_number_stats\n" );
	    break;
	case dis_star_data:
	    key = "DIS";
	    fprintf( fp, "save_generic_distance_number_stats\n" );
	    break;
	default:
            FATAL("Unknow case; code error ?");
	    exit(99); /* Will never be executed */
    }

    fprintf( fp, "   _Saveframe_category   distance_constraint_number_stats\n\n" );

    fprintf( fp, "   loop_\n" );
    fprintf( fp, "       _Constraint_list_label\n" );
    fprintf( fp, "       $distance_constraint_%s_%s\n", key, Mol_label );
    fprintf( fp, "\n   stop_\n\n" );
    
    if ( qts == noe_star_data )
    {
	fprintf( fp, "   _NOE_total_number                %6d\n", dist_count );
    
	strcpy( msg_string, "could not find index for contact range type" );
    
	fprintf( fp, "   _NOE_intraresidue_total_number   " );
	ix = ContactRangeIndex( intra_rng );
	if ( ix < 0 )
	{
	    ERROR( msg_string );
	    fprintf( fp, ".\n" ); 
	}
	else
	    fprintf( fp, "%6d\n", cont_typ[ix].tot_count ); 
    
	fprintf( fp, "   _NOE_sequential_total_number     " );
	ix = ContactRangeIndex( seq_rng );
	if ( ix < 0 )
	{
	    ERROR( msg_string );
	    fprintf( fp, ".\n" ); 
	}
	else
	    fprintf( fp, "%6d\n", cont_typ[ix].tot_count ); 
    
	fprintf( fp, "   _NOE_medium_range_total_number   " );
	ix = ContactRangeIndex( medium_rng );
	if ( ix < 0 )
	{
	    ERROR( msg_string );
	    fprintf( fp, ".\n" ); 
	}
	else
	    fprintf( fp, "%6d\n", cont_typ[ix].tot_count ); 
    
	fprintf( fp, "   _NOE_long_range_total_number     " );
	ix = ContactRangeIndex( long_rng );
	if ( ix < 0 )
	{
	    ERROR( msg_string );
	    fprintf( fp, ".\n" ); 
	}
	else
	    fprintf( fp, "%6d\n", cont_typ[ix].tot_count ); 
	
	ix = ContactRangeIndex( inter_rng );
	if ( ix < 0 )
	{
	    ERROR( msg_string );
	}
        else
/*	else if ( cont_typ[ix].tot_count != 0 )  always print this data */
	{
	    fprintf( fp, "   _NOE_intermolecular_total_number " );
	    fprintf( fp, "%6d\n", cont_typ[ix].tot_count ); 
	}
    }

    else if ( qts == hb_star_data )
    {
	fprintf( fp, "   _Constraints_hydrogen_bonds_constrained_total_number   @\n" );
	fprintf( fp, "   _Constraints_defining_hydrogen_bonds_total_number      %6d\n", dist_count );
    }

    else if ( qts == ss_star_data )
    {
	fprintf( fp, "   _Constraints_disulfide_bonds_constrained_total_number   @\n" );
	fprintf( fp, "   _Constraints_defining_disulfide_bonds_total_number      %6d\n", dist_count );
    }
    
    else
    {
	fprintf( fp, "   _Constraints_distance_assumed_total_number   %6d\n", dist_count );
    }

    fprintf( fp, "\nsave_\n\n" );
}

/****************  NewDrstcont  ****************/
/*
**  Allocate and initialize a Drstcont array.
*/

Drstcont *NewDrstcont( int count )
{
    Drstcont   *new_variable, *new_array;

    new_array = CREATE_NEW( Drstcont, count );
    if ( !new_array )
	NoMemErr( "NewDrstcont" );

    new_variable = new_array;
    while ( count-- )
    {
	new_variable->upp_count = 0;
	new_variable->low_count = 0;
	new_variable->uplo_count = 0;
	new_variable->tot_count = 0;
	new_variable++;
    }
    return( new_array );
}

/****************  UpdDrstcont  ****************/
/*
**  Update a Drstcont struct.
*/

void UpdDrstcont( Drstcont *cont, Distrst *rst )
{
    char *fct_name = "UpdDrstcont";

    switch ( rst->typ_bnd )
    {
	case 'B':
	case 'b':
	    cont->uplo_count++;
	    break;
	case 'u':
	    cont->upp_count++;
	    break;
	case 'l':
	    cont->low_count++;
	    break;
	default:
	    FATAL( "program/data error" );
    }
    cont->tot_count++;
}

/****************  WriteDrstcont  ****************/
/*
**  Write a Drstcont struct.
*/

void WriteDrstcont( FILE *fp, Drstcont *cont )
{
    fprintf( fp, " %6i %6i %6i %6i\n", 
             cont->low_count, cont->upp_count, cont->uplo_count, 
             cont->tot_count );
}

/****************  WriteHdrDrstcont  ****************/
/*
**  Write a header for subsequent calls to WriteDrstcont.
*/

void WriteHdrDrstcont( FILE *fp )
{
    fprintf( fp, " NLower NUpper NUp&Lo NTotal\n" );
}

/****************  WriteDstTyp  ****************/
/*
**  Write a distance type name.
*/

void WriteDstTyp( FILE *fp, int typ )
{
    fprintf( fp, "%4i ", typ );
}

/****************  UpdateDistInf  ****************/
/*
**  Update a Distinf list from old to new indices, using a map as
**  produced by RangeStrucset. (This is to update the Distinf list
**  to the new, shorter Strucset).
*/
 
void UpdateDistInf( Distinf *dist_inf, int dist_count, Mol *old_mol, Mol *new_mol, 
                    int *atom_map, int *residue_map )
{
    char      *fct_name = "UpdateDistInf";
    Residue   *new_residues = new_mol->residues;
    int       idx, skip;
    
    while ( dist_count-- )
    {

/* atom and residue A */

	skip = 0;
	idx = atom_map[dist_inf->idx_A];
	if ( idx == -1 )
	    skip++;
	dist_inf->idx_A = idx;
	idx = FindResidueIndexFromNum( dist_inf->res_A, old_mol );
	if ( idx < 0 )
	{
	    ERROR( "Index error - Code error" );
	    exit( 1 );
	}
	idx = residue_map[idx];
	if ( idx == -1 )
	{
	    skip++;
	    dist_inf->res_A = 0;
	}
	else
	    dist_inf->res_A = new_residues[idx].res_num_new;
	if ( skip == 1 )
	{
	    ERROR( "Inconsistent map info - Code error" );
	    exit( 1 );
	}

/* atom and residue B */

	skip = 0;
	idx = atom_map[dist_inf->idx_B];
	if ( idx == -1 )
	    skip++;
	dist_inf->idx_B = idx;
	idx = FindResidueIndexFromNum( dist_inf->res_B, old_mol );
	if ( idx < 0 )
	{
	    ERROR( "Index error - Code error" );
	    exit( 1 );
	}
	idx = residue_map[idx];
	if ( idx == -1 )
	{
	    skip++;
	    dist_inf->res_B = 0;
	}
	else
	    dist_inf->res_B = new_residues[idx].res_num_new;
	if ( skip == 1 )
	{
	    ERROR( "Inconsistent map info - Code error" );
	    exit( 1 );
	}

/* end of loop */

	dist_inf++;
    }
}

/****************  RangeDistLists  ****************/
/*
**  Apply a Range to a {Distinf, Distinfl, Distrst, Distval} list combination, 
**  i.e. select residues and/or models. New lists are created, and pointers
**  are returned.
**  (The difference with UpdateDistInf is that the indices remain the same, 
**  but certain restraints/models will be removed from the lists.)
**
**  If a related Distinfl list exists it is processed as well. It should
**  correspond to the normal list! (cf. "special case" processing in 
**  ProcessDistinfFile; also see main program.)
**
**  NOTE: Models are selected on the basis of struc_num, not fil_struc_num.
**  (The <models> array is currently allocated and filled by ProcessDistvalFile. 
**  A less ad hoc solution needs to consider the RangeStrucset routine.)
**  The Range is a global variable, used in routines InModelRange and 
**  InResidueRange.
*/

void RangeDistLists( DistdataList *inf_list, Distinf *dist_inf, Distrst **dist_rst,
                     DistdataList *val_list, int *dist_count, int *tot_dist_count,
		     int *a_count, Model **models, 
		     int *new_struc_count, const int struc_count, Mol *molptr )
{
    char     *fct_name = "RangeDistLists";
    int      d_count, s_count, old_count, new_count, idx_A, idx_B, how_sel;
    int      a_cnt = 0, tot_n_count = 0, count;
    Distinf  *new_inf, *n_inf, *d_inf;
    Distinfl *n_infl, *d_infl;
    Distinfl *new_infl, *dist_infl;
    Distrst  *new_rst, *n_rst, *d_rst;
    Distval  *dist_val, *new_val, *n_val, *n, *d_val, *v;
    Distvall *dist_vall, *new_vall, *n_vall, *nl, *d_vall, *vl;
    Model    *new_models, *n_mdl, *mdl;
    Boolean  *keep, *kp, select_A, select_B;
    Residue  *residues = molptr->residues;

    DistdataDesc   *inf_elem, *new_inf_elem, *new_val_elem;
    DistdataTyp    val_typ;

/* find the inf-DistdataList element */

    count = inf_list->count;
    inf_elem = inf_list->fst_desc;
    while ( count-- )
    {
	if ( inf_elem->data != dist_inf )
	    inf_elem = inf_elem->next;
	else
	    break;
    }
    
    dist_infl = NULL;
    if ( inf_elem == NULL )
	FATAL( "Could not find DistdataList element - code error" );
    if ( inf_elem->type != dd_inf_normal )
	FATAL( "Not normal type - code error" );
    if ( inf_elem->related )       /* we have a related (extended) list */
    {
	if ( inf_elem->related->type != dd_inf_extended ) {
	    FATAL( "Not extended type - code error" );
        } else {
	    dist_infl = (Distinfl *) (inf_elem->related->data);
        }
    }

/* set distval and determine its type */
/* NOTE: as with the inf list we are only looking at the last entry into
** the DistdataList. Later we may want to process multiple entries.
*/
    val_typ  = val_list->lst_desc->type;
    if ( val_typ == dd_val_normal )
    { 
	dist_val  = (Distval *) (val_list->lst_desc->data);
	dist_vall = NULL;
    }
    else if ( val_typ == dd_val_extended )
    {
	dist_val  = NULL;
	dist_vall = (Distvall *) (val_list->lst_desc->data);
    }
    else {
	FATAL( "Type not normal and not extended - code error" );
        exit(99); /* Never executed but for compiler warnings */
    }

/* count the number of models to be kept, and store their numbers */

    new_count = 0;
    new_models = CREATE_NEW( Model, struc_count );
    if ( !new_models )
	NOMEM( "new_models" );
    n_mdl = new_models;

    mdl = *models;
    s_count = struc_count;
    while ( s_count-- )
    {
	if ( InModelRange( mdl->struc_num ) )   /* use seq. not original number */
	{
	    new_count++;
	    n_mdl++->struc_num = mdl->struc_num;
	}
	mdl++;
    }
    
    if ( terseness < 6 ) {
        printf( " * Selected %i models\n", new_count );
    }
    if ( !new_count )
    {
	puts( " * Error" );
	exit( 1 );
    }
    
/* free the old list, and save the addres of the new one */

    free( *models );
    *models = new_models;
    *new_struc_count = new_count;
    
/* how to select restraints? */

    how_sel = GetOptionSubItem( opt_rest_sel );
    switch ( how_sel )
    {
	case restsel_both:
            if ( terseness < 6 ) {
	        puts( " * Selecting restraints for which BOTH residues are in allowed range" );
            }
	    break;
	case restsel_one:
            if ( terseness < 6 ) {
	        puts( " * Selecting restraints for which AT LEAST ONE residues is in allowed range" );
            }
	    break;
	default:
	    FATAL( "Illegal opt_rest_sel value - program error" );
    }

/* count the number of restraints to be kept */

    new_count = 0;
    old_count = *dist_count;
    keep = CREATE_NEW( Boolean, old_count );
    if ( !keep )
	NOMEM( "keep" );
    
    d_inf = dist_inf;
    kp    = keep;
    d_count = old_count;
    while ( d_count-- )
    {
	idx_A = FindResidueIndexFromNum( d_inf->res_A, molptr );
	idx_B = FindResidueIndexFromNum( d_inf->res_B, molptr );
	select_A = idx_A < 0 ? false : InResidueRange( residues+idx_A );
	select_B = idx_B < 0 ? false : InResidueRange( residues+idx_B );
	if (  how_sel == restsel_both ? ( select_A && select_B ) : ( select_A || select_B ) )
	{
	    new_count++;
	    *kp = true;
/*
            printf( "%i %i %i %i %i\n", d_inf->idx_A, d_inf->idx_B,
                    d_inf->res_A, d_inf->res_B, d_inf->dtype );
*/
	}
	else
	    *kp = false;
	d_inf++;
	kp++;
    }
    
    if ( terseness < 6 ) {
        printf( " * Selected %i restraints\n", new_count );
    }
    if ( !new_count )
    {
	puts( " * Error" );
	exit( 1 );
    }

/* allocate new lists */

    new_inf = CREATE_NEW( Distinf, new_count );
    if ( dist_infl )
	new_infl = CREATE_NEW( Distinfl, new_count );
    else
        new_infl = NULL; /* Will not be used */
    
    new_rst = CREATE_NEW( Distrst, new_count );
    
    /* use explicite initializations against compiler warnings */
    if ( val_typ == dd_val_normal ) 
    {
	new_val    = CREATE_NEW( Distval, new_count * struc_count );
        new_vall   = NULL;
	n_val      = new_val;
	n_vall     = NULL;
	n          = new_val;
        nl         = NULL;
    }
    else
    {
        new_val    = NULL;
	new_vall   = CREATE_NEW( Distvall, new_count * struc_count );
	n_val      = NULL;
	n_vall     = new_vall;
        n          = NULL;
	nl         = new_vall;
    }
    if ( !new_val && !new_vall )
	NOMEM( "new_inf/new_rst/new_val" );
    
    n_inf = new_inf;
    n_rst = new_rst;


/* copy the info to the new lists */

    d_inf = dist_inf;
    kp    = keep;
    d_rst = *dist_rst;
    d_count = old_count;
    if ( val_typ == dd_val_normal )
    {
	d_val = dist_val;
	v     = dist_val;
	while ( d_count-- )
	{
	    if ( *kp )
	    {
		*n_rst++ = *d_rst;
		*n_inf++ = *d_inf;
		mdl = *models;
		for ( s_count = 1; s_count <= struc_count; s_count++ )
		{
		    if ( mdl->struc_num == s_count )
		    {
			mdl++;
			*n = *v; 
			n += new_count;
		    }
		    v += old_count;
		}
		n = ++n_val;
	    }
	    v = ++d_val;
	    ++d_rst;
	    ++d_inf;
	    ++kp;
	}
    }
    else
    {
	d_vall = dist_vall;
	vl     = dist_vall;
	while ( d_count-- )
	{
	    if ( *kp )
	    {
		*n_rst++ = *d_rst;
		*n_inf++ = *d_inf;
		tot_n_count += vl->count;
		mdl = *models;
		for ( s_count = 1; s_count <= struc_count; s_count++ )
		{
		    if ( mdl->struc_num == s_count )
		    {
			mdl++;
			*nl = *vl; 
			nl += new_count;
		    }
		    vl += old_count;
		}
		nl = ++n_vall;
	    }
	    vl = ++d_vall;
	    ++d_rst;
	    ++d_inf;
	    ++kp;
	}
    }


/* and now for the new Distinfl list */

    if ( dist_infl )
    {
	n_infl = new_infl;
	d_infl = dist_infl;
	kp     = keep;
	d_count = old_count;
	while ( d_count-- )
	{
	    if ( *kp )
	    {
		a_cnt += ( d_infl->count_A + d_infl->count_B );
		*n_infl++ = *d_infl;
	    }
	    ++d_infl;
	    ++kp;
	}
    }

/* free the old lists, and save the addresses of the new ones */
    
    free( keep );
    free( *dist_rst );
    *dist_rst = new_rst;
    *dist_count = new_count;
    *tot_dist_count = tot_n_count;
    *a_count = a_cnt;
    

    DestroyDistdataListElem( inf_list, dist_inf, true );
    new_inf_elem = NewDistdataListElem( inf_list, NULL );
    if ( !dist_infl )
    {
	new_inf_elem->data = (void *)new_inf;
	new_inf_elem->type = dd_inf_normal;
    }
    else          /* keep the same order as in ProcessDistinfFile */
    {
	new_inf_elem->data = (void *)new_infl;
	new_inf_elem->type = dd_inf_extended;
        new_inf_elem = NewDistdataListElem( inf_list, inf_list->lst_desc );
	new_inf_elem->data = (void *)new_inf;
	new_inf_elem->type = dd_inf_normal;
    }
    DestroyDistdataListElem( val_list, dist_val, true );
    new_val_elem = NewDistdataListElem( val_list, NULL );
    if ( val_typ == dd_val_normal )
	new_val_elem->data = new_val;
    else
	new_val_elem->data = new_vall;
    new_val_elem->type = val_typ;
}

/****************  ClassifyRestraints  ****************/
/*
**  Classify restraints in a Distinf list according to residue range.
**  Cf. code in SortRestraints.
*/

void ClassifyRestraints( Distinf *dist_inf, int dist_count, Atom *atoms )
{
    char           *fct_name = "ClassifyRestraints";
    ContactRange   contact_range;
    int            diff_resnum, insert_present = false, nerr = 0;

    while ( dist_count-- )
    {
	if ( atoms[dist_inf->idx_A].res_insert != ' ' || atoms[dist_inf->idx_B].res_insert != ' ' )
	    insert_present = true;

	diff_resnum = atoms[dist_inf->idx_A].res_num_orig - atoms[dist_inf->idx_B].res_num_orig ;
	
	if ( diff_resnum < 0 )
	    diff_resnum *= -1; 
	   	
	contact_range = undef_rng;
	if ( diff_resnum == 0 )
	    contact_range = intra_rng;
	if ( diff_resnum == 1 )
	    contact_range = seq_rng;
	if ( ( diff_resnum > 1 ) && ( diff_resnum <= max_MEDIUMRANGE ) )
	    contact_range = medium_rng;
	if ( diff_resnum > max_MEDIUMRANGE )
	    contact_range = long_rng;
	if ( NOTEQID( atoms[dist_inf->idx_A].chain, atoms[dist_inf->idx_B].chain ) )
	    contact_range = inter_rng;
	
	if ( contact_range == undef_rng )
	    nerr++;
	dist_inf->dtype = (int)contact_range;
	dist_inf++;
    }
    if ( insert_present )
	ERROR( "classification partially incorrect due to presence of residue insertion code(s)" );
    if ( nerr )
    {
	sprintf( msg_string, "%i restraints could not be classified according to range", nerr );
	ERROR( msg_string );
    }
}

/****************  ShowDistl  ****************/
/*
** Print all data from a set of distinf/distval/distrst arrays (long version)
*/

void ShowDistl( Distinfl *dist_inf, Distvall *dist_val, Distrst *dist_rst,
                const int dist_count, Mol *molptr, const int struc_count )
{
    FILE       *fp = stdout;
    Distinfl   *infptr;
    Distvall   *valptr, *vptr;
    Distrst    *rstptr;
    Distidx    *A, *B;
    Atom       *atoms = molptr->atoms;
    int        count_A, count_B, count, countStr, num;
   
    infptr = dist_inf;
    valptr = dist_val;
    rstptr = dist_rst;

    count = dist_count;
    while ( count-- )
    {
	if ( valptr->count != ( infptr->count_A * infptr->count_B ) )
	{
	    puts( " * inconsistent count" );
	    exit( 1 );
	}
        num = 0;
	A       = infptr->A;
	count_A = infptr->count_A;
	while ( count_A-- )            /* loop over atom A */
	{
	    B       = infptr->B;
	    count_B = infptr->count_B;
	    while ( count_B-- )        /* loop over atom B */
	    {
		WriteFAtom( fp, &atoms[ A->idx ], use_ext );
		fputc( ' ', fp );
		WriteFAtom( fp, &atoms[ B->idx ], use_ext );
		fprintf( fp, " <%6.2f>", rstptr->up_bnd );

		vptr   = valptr;       /* Distvall belonging to this restraint */
                                       /* first model */
		countStr = struc_count;
		while ( countStr-- )   /* loop over models */
		{

/* now we have to compute the proper location in the values array     */
/* <num> is also equal to [i*count_B + j], where i and j are the indices */
/* in A and B:  0<=i<count_A and 0<=j<count_B */

		    fprintf( fp, " %6.2f", vptr->values[num] );
		    vptr += dist_count;
                                       /* jump forward in Distvall array */
                                       /* to next structure */
		}
		fputc( '\n', fp );
		num++;
		B++;
	    }
	    A++;
	}
	fputc( '\n', fp );
	infptr++;
	valptr++;        /* now increment to the next Distvall element */
        rstptr++;
    }
}

